

In [1]:

    
import camb
pars = camb.CAMBparams()
from scipy.special import jn, jn_zeros
from camb import model, initialpower
from scipy.interpolate import interp1d
from hankel_transform import *
from astropy.cosmology import Planck13 #use Planck15 if you can
import astropy.units as u
#rc('text', usetex=False)



In [2]:

    
from matplotlib.pyplot import *
%matplotlib inline



In [3]:

    
cosmo=Planck13.clone()#we want h=1 for some comparisons



In [4]:

    
#Set up a new set of parameters for CAMB
kmax=55#30
kmin=.8e-3
k_smooth=1
nk=5000
rmin=.055
rmax=110
non_linear=1

pars = camb.CAMBparams()
#This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency
pars.set_cosmology(H0=67.27, ombh2=0.022250, omch2=0.119800, mnu=0.06, omk=0, tau=0.06)
pars.InitPower.set_params(ns=0.965, r=0,As =2.14e-09)
pars.set_for_lmax(2500, lens_potential_accuracy=0)
zb=[.27]#[0.0]
pars.set_matter_power(redshifts=zb,kmax=kmax);
if non_linear==1:
    pars.NonLinear = model.NonLinear_both
else:
    pars.NonLinear = model.NonLinear_none
results = camb.get_results(pars)
kh, z, pk =results.get_matter_power_spectrum(minkh=kmin, maxkh=kmax, npoints =nk)



In [5]:

    
#Setting up the Hankel Transform
#This part is slower. But only needs to be run once. 
#If you only need wgg, set j_nu=[0]. For wg+ (or \Delta\Sigma) use j_nu=[2]
%time HT=hankel_transform(rmin=rmin,rmax=rmax,kmax=kmax,j_nu=[0],n_zeros=28000,kmin=kmin)
#HT=hankel_transform(rmin=1,rmax=rmax,kmax=1,j_nu=[0,2],n_zeros=2800,kmin=1.e-2)#quick test... inaccurate









    



('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 29000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 30000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 31000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 32000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 33000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 34000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 35000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 36000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 37000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 38000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 39000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 40000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 41000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 42000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 43000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 44000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 45000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 46000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 47000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 48000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 49000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 50000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 51000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 52000)
('j-nu=', 0, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 53000)
('nr:', 1927)
CPU times: user 53.3 s, sys: 715 ms, total: 54 s
Wall time: 54.3 s



In [6]:

    
pk_taper=HT.taper(k=kh,pk=pk[0],large_k_lower=5,large_k_upper=kmax,low_k_lower=kmin,
                     low_k_upper=kmin*2)
#need atleast k=10 and k=1.e-3 to get decent wgg for 1-100 Mpc/h range.

r_gg,wgg=HT.projected_correlation(k_pk=kh,pk=pk[0]*b_g**2,j_nu=0) %time r_gg,wgg_taper=HT.projected_correlation(k_pk=kh,pk=pk_taper*b_g**2,j_nu=0)

b_g=1

#Internal test case wgg_test_fact=1./np.pi**2 wgg_test=np.genfromtxt('../data/planck15_wgg_rsd0_kz0_s0_k10_z0.27.dat',names=('rp','wgg')) wgp_test=np.genfromtxt('../data/planck15_wgp_rsd0_kz0_s0_k10_z0.27.dat',names=('rp','wgp'))

plot(kh,pk[0]) plot(kh,pk_taper) xscale('log') yscale('log') xlim(1.e-4,30) ylim(1,5.e4)

plot(r_gg,r_gg*wgg,label='No tapering') plot(r_gg,r_gg*wgg_taper,label='tapering') #plot(wgg_test['rp'],wgg_test['rp']*wgg_test['wgg']*wgg_test_fact) xscale('log') #yscale('log') #ylim(1.e-2,1.e3) xlim(1,110) xlabel('$r_p$ [Mpc/h]') ylabel('$w_{gg}$ [Mpc/h]') ylim(0,72) legend()

plot(r_gg,wgg-wgg_taper,'r') xscale('log') xlabel('$r_p$ [Mpc/h]')

w_intp=interp1d(r,w,bounds_error=False,fill_value=np.nan) plot(wgg_test['rp'],w_intp(wgg_test['rp'])/(wgg_test['wgg']*wgg_test_fact),'r--') w_intp=interp1d(r,w2,bounds_error=False,fill_value=np.nan) plot(wgg_test['rp'],w_intp(wgg_test['rp'])/(wgg_test['wgg']*wgg_test_fact)) xscale('log') ylim(.98,1.02)

g-g Lensing covariance



In [7]:

    
h = 0.704



In [8]:

    
pars = camb.CAMBparams()
#This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency

# replaced with DarkSky values
pars.set_cosmology(H0=70.4, ombh2=0.0231566, omch2=0.12293, mnu=0.06, omk=0, tau=0.06)
pars.InitPower.set_params(ns=0.967, r=0,As =2.18e-09)
pars.set_for_lmax(2500, lens_potential_accuracy=0)
zb=[0.0]#[0.0]#
pars.set_matter_power(redshifts=zb,kmax=kmax);
if non_linear==1:
    pars.NonLinear = model.NonLinear_both
else:
    pars.NonLinear = model.NonLinear_none
results = camb.get_results(pars)
kh, z, pk =results.get_matter_power_spectrum(minkh=kmin, maxkh=kmax, npoints =nk)



In [9]:

    
HT=hankel_transform(rmin=rmin,rmax=rmax,kmax=kmax,j_nu=[2],n_zeros=60000,kmin=kmin)









    



('changed kmax to', 93.37495094255787, ' to cover rmin')
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 61000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 62000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 63000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 64000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 65000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 66000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 67000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 68000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 69000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 70000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 71000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 72000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 73000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 74000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 75000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 76000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 77000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 78000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 79000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 80000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 81000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 82000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 83000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 84000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 85000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 86000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 87000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 88000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 89000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 90000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 91000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 92000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 93000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 94000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 95000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 96000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 97000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 98000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 99000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 100000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 101000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 102000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 103000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 104000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 105000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 106000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 107000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 108000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 109000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 110000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 111000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 112000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 113000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 114000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 115000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 116000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 117000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 118000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 119000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 120000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 121000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 122000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 123000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 124000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 125000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 126000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 127000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 128000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 129000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 130000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 131000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 132000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 133000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 134000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 135000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 136000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 137000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 138000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 139000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 140000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 141000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 142000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 143000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 144000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 145000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 146000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 147000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 148000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 149000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 150000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 151000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 152000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 153000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 154000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 155000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 156000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 157000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 158000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 159000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 160000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 161000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 162000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 163000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 164000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 165000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 166000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 167000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 168000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 169000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 170000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 171000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 172000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 173000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 174000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 175000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 176000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 177000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 178000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 179000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 180000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 181000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 182000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 183000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 184000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 185000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 186000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 187000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 188000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 189000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 190000)
('j-nu=', 2, ' not enough zeros to cover kmin, increasing by ', 1000, ' to', 191000)
('nr:', 3269)



In [10]:

    
sigma_crit=4.7e3



In [11]:

    
sigma_e=.36
sigma_gamma=sigma_e/1.7
n_s=8
shape_noise=sigma_crit**2*sigma_gamma**2/n_s#*cosmo.H_z(z=0.27)/cosmo.c



In [12]:

    
rho=0.04180594698596614 #units of Msun/pc^2/mpc, using h=1



In [13]:

    
Dchi2=900 #integral over product of window functions
p_kappa_kappa=pk[0]*(rho)**2*Dchi2



In [14]:

    
n_g=5e-4#3.e-4
g_shot_noise=1./n_g
b_g=1.838
p_g=b_g**2*pk[0]
#p_gk=b_g*pk[0]*rho



In [15]:

    
r_bins=np.logspace(-1,1.6,19)#*h



In [16]:

    
area=10000
area_comoving=area*(np.pi/180)**2*cosmo.comoving_distance(z=.27)**2
L_W=500
#vol=(area_comoving*L_W).value
vol=(1000)**3#vol.value # don't actually need the h



In [17]:

    
(shape_noise*g_shot_noise)/(vol*h**4)









    Out[17]:





1.0082128926958325

vol/vol_old.value

from pearce.mocks.kittens import DarkSky cat = DarkSky(264, system = 'sherlock')

#real_p_g = np.load('/home/users/swmclau2/Git/pearce/bin/covmat/galaxy_pk.npy') real_p_g = np.load('/home/users/swmclau2/Git/pearce/notebooks/p_g.npy')

real_p_g = real_p_g[1:]



In [18]:

    
taper_kw=dict({'large_k_lower':10,'large_k_upper':kmax,'low_k_lower':kmin,'low_k_upper':kmin*1.2})



In [19]:

    
r,cov_ggkk_shape=HT.projected_covariance(k_pk=kh,pk1=p_g+g_shot_noise,pk2=np.zeros_like(p_kappa_kappa)+shape_noise,j_nu=2,taper=True,**taper_kw)
r,cov_ggkk=HT.projected_covariance(k_pk=kh,pk1=p_g+g_shot_noise,pk2=p_kappa_kappa,j_nu=2,taper=True,**taper_kw)

r_re,cov_ggkk_shape_re=HT.bin_cov(r=r/h,cov=cov_ggkk_shape,r_bins=r_bins)

r_re,cov_ggkk_re=HT.bin_cov(r=r/h,cov=cov_ggkk,r_bins=r_bins)
corr=HT.corr_matrix(cov=cov_ggkk_shape_re)

cov_ggkk_re/=vol # Don't forget!
cov_ggkk_shape_re/=vol



In [20]:

    
print np.sqrt(np.diag(cov_ggkk_shape_re))









    



[1.0263207  0.98494148 0.90722623 0.69274081 0.57677226 0.45103709
 0.32990732 0.23736836 0.17742547 0.13398334 0.10010477 0.07595037
 0.05802125 0.04515957 0.03575584 0.02910926 0.02429026 0.02063339]

# I've corrected this so its no longer necessary cov_ggkk_shape_re/=(h**4) # convert to the units I'm in cov_ggkk_re/=(h**4)



In [21]:

    
print r_re*np.sqrt(np.diag(cov_ggkk_shape_re))









    



[0.12120075 0.1622098  0.20836593 0.22188417 0.2576347  0.28096789
 0.28660307 0.28757834 0.29977384 0.31569886 0.32894359 0.34804948
 0.37080248 0.40248509 0.44441775 0.50456842 0.5871721  0.69558191]

plot(kh, pk.squeeze()) plot(real_p_g[:,0], real_p_g[:,1]) loglog();

np.sqrt(real_p_g[np.argmin(np.abs(real_p_g[:,0] - 0.1) ),1]/pk.squeeze()[np.argmin(np.abs(kh - 0.1) )])

r,cov_ggkk=HT.projected_covariance(k_pk=real_p_g[:,0],pk1=real_p_g[:,1]+g_shot_noise,pk2=np.zeros_like(real_p_g[:,0]) + shape_noise,j_nu=2,taper=True,**taper_kw) r_re,cov_ggkk_re=HT.bin_cov(r=r,cov=cov_ggkk,r_bins=r_bins) corr=HT.corr_matrix(cov=cov_ggkk_re)

np.sqrt(np.diag(cov_ggkk_re))



In [22]:

    
sim_covmat = np.load('/home/users/swmclau2/Git/pearce/bin/covmat/wp_ds_sim_covmat.npy')

plot(r_re, np.sqrt(np.diag(cov_ggkk_re)) ) plot(r_re, np.sqrt(np.diag(cov_ggkk_shape_re)), color ='r' ) plot(r_re, np.sqrt(np.diag(sim_covmat[18:][:,18:])), color = 'g') ylabel('Shape Noise Error') xlabel('r [Mpc]') loglog();



In [24]:

    
plot(r_re, np.ones_like(np.sqrt(np.diag(cov_ggkk_re)) ) )
#plot(r_re, np.sqrt(np.diag(cov_ggkk_shape_re)), color ='r' )
#plot(r_re, np.sqrt(np.diag(sim_covmat[18:][:,18:]))/np.sqrt(np.diag(cov_ggkk_re)), color = 'g')
ylabel('Shape Noise Error')
xlabel('r [Mpc]')
loglog();









    



/share/software/user/open/py-scipystack/1.0_py27/lib/python2.7/site-packages/matplotlib/axes/_base.py:3179: UserWarning: Attempting to set identical bottom==top results
in singular transformations; automatically expanding.
bottom=1.0, top=1.0
  'bottom=%s, top=%s') % (bottom, top))



In [26]:

    
#np.sqrt(np.diag(sim_covmat[18:][:,18:]))/np.sqrt(np.diag(cov_ggkk_re))



In [27]:

    
#plot(r_re, np.sqrt(np.diag(cov_ggkk_re)) )
plot(r_re, r_re*np.sqrt(np.diag(cov_ggkk_shape_re)), color ='r' )

ylabel('Shape Noise Error')
xlabel('r [Mpc]')
loglog();



In [28]:

    
print np.sqrt(np.diag(corr))









    



[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]



In [29]:

    
np.sqrt(np.diag(cov_ggkk_shape_re))









    Out[29]:





array([1.0263207 , 0.98494148, 0.90722623, 0.69274081, 0.57677226,
       0.45103709, 0.32990732, 0.23736836, 0.17742547, 0.13398334,
       0.10010477, 0.07595037, 0.05802125, 0.04515957, 0.03575584,
       0.02910926, 0.02429026, 0.02063339])



In [30]:

    
imshow(corr)









    Out[30]:





<matplotlib.image.AxesImage at 0x7faaf452d350>



In [31]:

    
np.save('shape_noise.npy', cov_ggkk_shape_re)



In [32]:

    
np.min(r_re), np.min(r)









    Out[32]:





(0.11809247432396976, 0.055)



In [33]:

    
from matplotlib.colors import LogNorm



In [34]:

    
pcolor(r_re,r_re,corr,norm=LogNorm())
colorbar()
xscale('log')
yscale('log')
xlim(HT.rmin,HT.rmax)
ylim(HT.rmin,HT.rmax)









    Out[34]:





(0.055, 110)



In [35]:

    
r,cov_gkgk=HT.projected_covariance(k_pk=kh,pk1=p_gk,pk2=p_gk,
                               kmax=100,rmin=.8,j_nu=2,rmax=110,n_zeros=3500)#return_Jrr=True,Jrr=Jrr
cov_gkgk*=Dchi2
r_re,cov_gkgk_re=HT.bin_cov(r=r,cov=cov_gkgk,r_bins=r_bins)
corr=HT.corr_matrix(cov=cov_gkgk_re)









    




NameErrorTraceback (most recent call last)
<ipython-input-35-98df143096fb> in <module>()
----> 1 r,cov_gkgk=HT.projected_covariance(k_pk=kh,pk1=p_gk,pk2=p_gk,
      2                                kmax=100,rmin=.8,j_nu=2,rmax=110,n_zeros=3500)#return_Jrr=True,Jrr=Jrr
      3 cov_gkgk*=Dchi2
      4 r_re,cov_gkgk_re=HT.bin_cov(r=r,cov=cov_gkgk,r_bins=r_bins)
      5 corr=HT.corr_matrix(cov=cov_gkgk_re)

NameError: name 'p_gk' is not defined



In [ ]:

    
cov_final=(cov_ggkk_re+cov_gkgk_re)/vol
corr=HT.corr_matrix(cov=cov_final)
errors=HT.diagonal_err(cov=cov_final)



In [ ]:

    
pcolor(r_re,r_re,corr,norm=LogNorm())
colorbar()
xscale('log')
yscale('log')
xlim(HT.rmin,HT.rmax)
ylim(HT.rmin,HT.rmax)

error_qpm=np.genfromtxt('../data/qpm_DS_error_th.dat',names=True)



In [ ]:

    
plot(r_re,r_re*errors,label='new')
#plot(error_qpm['rp'],error_qpm['rp']*error_qpm['DS_err_gR'])
xscale('log')
yscale('log')
legend()
xlim(1,100)



In [ ]:



In [ ]: